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Abstract 

We present a study of the deconfinement phase transition of one-flavour QCD, using the multi- 
boson algorithm. The mass of the Wilson fermions relevant for this study is moderately large and 
the non-hermitian multiboson method is a superior simulation algorithm. Finite size scaling is 
studied on lattices of size 8 3 x 4, 12 3 x 4 and 16 3 x 4. The behaviours of the peak of the Polyakov 
loop susceptibility, the deconfinement ratio and the distribution of the norm of the Polyakov loop 
are all characteristic of a first-order phase transition for heavy quarks. As the quark mass de- 
creases, the first-order transition gets weaker and turns into a crossover. To investigate finite size 
scaling on larger spatial lattices we use an effective action in the same universality class as QCD. 
This effective action is constructed by replacing the fermionic determinant with the Polyakov loop 
identified as the most relevant Z(3) symmetry breaking term. Higher-order effects are incorpo- 
rated in an effective Z(3)-breaking field, h, which couples to the Polyakov loop. Finite size scaling 
determines the value of h where the first order transition ends. Our analysis at the end - point, 
h ep , indicates that the effective model and thus QCD is consistent with the universality class of the 
three dimensional Ising model. Matching the field strength at the end point, h ep , to the k values 
used in the dynamical quark simulations we estimate the end point, n ep , of the first-order phase 
transition. We find K ep ~ 0.08 which corresponds to a quark mass of about 1.4 GeV . 

Keywords: Lattice QCD, Dynamical quarks, Deconfinement phase transition. 

PACS numbers: 11.15.Ha, 12.38.Ge, 12.38.Gc, 12.38.Mh, 64.60.Fr. 



1. Introduction 

Understanding the properties of Quantum Chromodynamics (QCD) under extreme conditions of high 
temperature and/or pressure is a challenging problem and was considered long ago. Polyakov [PJ 
and Susskind Q predicted that QCD will undergo a deconfinement phase transition from normal 
hadronic matter to quark-gluon plasma when the temperature is increased. Such a phase transition is 
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believed to have occurred, in the opposite direction, lfis after the Big Bang and its nature is therefore 
important in Astrophysics. Theoretical information on the deconfinement phase transition has also 
become important as planned ultra-relativistic experiments will soon start at RHIC, Brookhaven and 
later on at LHC, CERN. In these experiments the temperature reached will be of the order of 600 
MeV and, at this temperature, one is still dealing with a strongly interactive system. Thus Lattice 
QCD provides the most suitable non-perturbative approach to study such phenomena using directly 
the QCD Lagrangian. 

The zero flavour sector of the theory (quenched) has been studied extensively and it is established 
that there is a first - order deconfinement phase transition || with the critical temperature determined 
in the continuum limit [||, |j| ||, [?], ^. State of the art lattice calculations are now being done including 
pair creation. For two flavours (Nj = 2) there are simulations both with Wilson || and staggered 
fermions [||, whereas for greater number of flavours the simulations are usually performed with the 
standard Wilson action [10, O]. On the other hand, one - flavour QCD has been largely ignored (early 
exceptions are given in ref. |12{| ), perhaps because of algorithm difficulties, although it has interesting 
properties. The continuum Nf = 1 theory contains no pions and it is expected to have no chiral 
phase transition [13] with the results of a study of one flavour staggered fermion [14| implying a chiral 
phase transition coming as a surprise. The purpose of this work is to fill this gap by investigating the 
deconfinement phase transition of one - flavour QCD. 

A model that plays an important role in our understanding of the phase diagram in QCD is the 
three states Potts model in three dimensions. It has a Z(3) symmetry and the spontaneous break 
down of this symmetry is expected to drive the deconfinement phase transition like it does in quenched 
QCD. In the presence of an external field the Z(3) symmetry is explicitly broken just like the fermionic 
determinant breaks Z(3) symmetry in QCD. In fact it was shown by DeGrand and DeTar [15] that, 
at high temperature and heavy quark mass, QCD reduces to the 3 - states Potts model in an external 
field. From its phase diagram shown schematically in Fig. [I], we observe that the first order phase 
transition gets weaker as the strength of the external field h increases and ends at some critical value 
h ep with a second order transition. Beyond this point a crossover behaviour is seen. Whereas Fig. |l] 
gives us a good starting point for the expected qualitative behaviour for QCD with dynamical quarks, 
the quantitative question of existence of such an end-point and its location can only be answered after 
a detailed calculation. 

In the real world of two light and one heavy quarks the phase diagram depends crucially on the 
size of the light and heavy quark masses. Therefore it is usual to consider the phase diagram of 3- 
flavour QCD in the mass plane m s vs m u ^. The mass point (0,0) is the chiral limit with a first order 
transition. Pure gauge is the opposite limit of infinite quark masses where a first order transition 
is also established. As one moves away from the infinite mass (or hopping parameter k = 0) limit 
the first order transition is expected, like in the Potts model, to persist but to become weaker and 
eventually to disappear, as the quark mass is reduced. This presumed robustness is based on the 
expectation that an infinitesimal increase in k does not cause the finite gap in the thermodynamical 
observables to change discontinuously but smoothly. That a first order transition is in general robust 
is supported from our experience with spin systems where a first order transition remains first order for 
small magnetic fields. In this work we seek to determine the point on the line (rn u d = oo,m s ) where 
the first order transition ends as well as the universality class of the second order phase transition 
at this end point. It turns out that the relevant quark masses for this study are relatively large and 
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Figure 1: Phase diagram in the three states Potts model in three dimensions 

therefore the multiboson algorithm is a very efficient method for simulating such dynamical quarks. 

From the determination of the end point for one - flavour we can draw conclusions about the end 
point for two flavours i.e. where on the line (m u d, m s = oo) the first order transition ends. The quark 
mass where this is expected to occur is about hal the mass obtained in the one - flavour case because 
the quark effects that we observe are about half as large as the effects obtained for two flavours of the 
same mass. 

In order to determine quantitatively the phase diagram for QCD we perform simulations on lattices 
with size 8 3 x 4, 12 3 x 4 and 16 3 x 4. Performing a finite size scaling analysis we find indeed that 
there is a first order phase transition for heavy quarks which gets weaker as the quark mass decreases 
and becomes second order at a value of K ep ~ 0.08. For smaller quark masses a crossover behaviour 
is seen [jO]]. 

In order to determine the end point in the (/3, k) plane accurately one needs larger lattice sizes, 
which takes very long to simulate. However there is a different route that we can follow. Namely 
we investigate an effective model in the same universality class as QCD which is easier to simulate, 
determine the end point there and then match back to QCD. Such an effective model can be con- 
structed by taking the pure gauge action and adding the most relevant Z(3) breaking term which is 
a Polyakov loop. The strength of the Z(3) breaking field is adjusted to model higher order terms. 
Within this effective model we were able to perform a finite - size scaling analysis up to a lattice size 
of 24 3 and determine the critical strength of the Z(3) breaking field where the first order transition 
ends. We also investigated the lattice spacing dependence within the effective model, by performing 
a simulation for temporal extension Nf = 2 in addition to Nf = 4. Finite size scaling at the end - 
point, h ep , yields results that are consistent with the scaling behaviour seen in the three dimensional 
Ising model. We obtain h(n) by performing a best fit of the Polyakov loop histograms obtained in the 
effective model to those obtained in QCD. This non perturbative matching of h ep yields the critical 
value of K ep and of f3 ep . 

This paper is organized as follows: In section (2) we give details on the local bosonic algorithm 
which we used to simulate one dynamical Wilson fermion. In section (3) we describe the observables 
that we used to probe a change of phase and the nature of this phase transition. We give the results of 
this analysis in section (4). In section (5) we discuss the simulation and the finite size scaling analysis 



3 



of the effective model and connect it to full QCD. In section (6) we discuss the continuum limit and 
finally in section (7) we summarize and conclude. 



2. Local bosonic algorithm for one flavour 

The local bosonic algorithm was originally proposed by Liischer [17] as an alternative method to the 
widely used Hybrid Monte Carlo (HMC) algorithm to simulate dynamical quarks. The basic idea is 
the replacement of the fermionic determinant by a functional integral over n bosonic fields having a 
local action. If Q = 75.0/(1 + 8 k) where D is the fermionic Wilson matrix then for two degenerate 
flavours we have 

detQ 2 oc lim / TT #t#fce~S=i 4[«-^) 2 +* ( 2 .1) 

J k=l 

where Jz^ = /J>k + i> v k with z k being the roots of a polynomial of even degree n constructed so that 

lim P n (Q 2 ) = -1 . (2.2) 

n— >oo 



The error introduced by taking n finite can be eliminated by a global accept/reject Metropolis test [18]. 
The generalization to any number of flavours is made possible by finding a polynomial approximation 
to the fermionic matrix itself rather than to Q 2 ]18| , |19| ], This can be done by constructing a polynomial 
of even degree n defined in the complex plane with complex conjugate roots z k such that 

lim P n (z) = - (2.3) 

n^oc Z 

for any z in the domain containing the spectrum of D (not including the origin). Since the spectral 
radius of the hopping matrix M is bounded by 8 in the free case and even less in the interacting case, 
we are guaranteed that the spectrum of O = 1 — kM will remain in the complex right half-plane for 
the heavy to moderately heavy quarks that we simulate (k < 1/8). This is demonstrated in Fig. @ 
where the boundary of the spectrum of the Dirac matrix is estimated for k = 0.1 as in ref. [PQ] . 
Using the property O = 75-0^75 we obtain 

n/2 

detP n (D) = c n Yl det(D - z k )Uet(D - z k ) (2.4) 

k=l 



with c n an easily computed constant IS]. Instead of eq. ( |2.ip one now finds 

detD = lin^det- 1 [Tl /2 (D)T n/2 (D)] (2.5) 

n/2 

oc lim / n ^t# fc e"^i^ (D ^ fc)t(D " 2fc)0ft , 



rfi /2 

where T n / 2 (D) = Yl k=1 {D — z k ). The algorithm is made exact with a global Metropolis test as follows. 
The one-flavour determinant can be expressed as 



detO = Cdet" 1 



Tl /2 {D)T n/2 {D) (2.6) 
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Figure 2: Estimated boundary of the Dirac eigenvalue spectrum 
with the correction factor C given by 



C 



lim 

m— >oo 



det" 1 


Tl /2 (D)T m/2 (D) 


det" 1 


T ] n/2 {D)T n/2 {D) 



(2.7) 



Because the error of the polynomial approximation ( p.3[ ) decreases exponentially with the degree of 
the polynomial, it is not necessary to actually take the limit m — > oo above. We observed that 
taking m > 3n, where n is adjusted for sufficient Metropolis acceptance below, is enough to make the 
systematic approximation error completely negligible compared to the statistical one. We implemented 



the correction factor C by a noisy Metropolis test, along the lines suggested in ref.|]18|. The gauge 
and boson fields (U, 4>) are updated by a sequence of local Monte Carlo steps forming a trajectory 
(U,4>) — ► (U 1 , ft). This procedure satisfies detailed balance with respect to the approximate action 



ra/2 



-"approx 



S g [U] + J2\( D - z k)<Pk 



k=l 



At the end of each trajectory, one generates a field r\ with probability 

P HB ( V ) = Re-W 2 

where X = T m j 2 (D)T~h(D) and R is a normalization constant. The new configuration (U', 
accepted with probability 



(2- 



IS 



pA 



mm 



1, 



-\X' V \' 



(2.10) 



The correction term is thus estimated by only one r\ field. Detailed balance with respect to the 
desired action det -1 TL /2 {p)T m i 2 (D) e~ Sa ^ is satisfied after averaging the probability density for 



m/2\ 

the Metropolis test, P HB (rj)P ( j, 



(£/,0)-(£/' ,</>') 



a, over the r\ field. 
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The number of bosonic fields n is chosen so that the correction term leads to an acceptance rate 
of about 2/3 and m is taken at least three times n. 

For the local updating of the gauge and boson fields we used standard heatbath and over-relaxation 



algorithms as described in [18|. A trajectory is a symmetric combination of (21 + 1) over-relaxation 
steps applied alternatively to the gauge and boson fields, preceded and followed by a heatbath on the 
bosons. Ergodicity for the gauge fields is maintained due to their coupling to the bosonic fields. The 
roots Zk are distributed on a circle centered at (1,0). We implemented even - odd preconditioning to 
lower the number of bosonic fields needed for a given accuracy. To efficiently equilibrate the system, 
we start from thermalized quenched gauge configurations and initialize the boson fields by generating 
a gaussian random vector x and setting 

n 

<Pk^c n D]l(D-z k )x ■ (2.11) 

3. Signals for the phase transition or lack thereof 

QCD with infinitely heavy quarks, i.e. quenched, undergoes a first-order deconfinement transition 



corresponding to the spontaneous breaking of the Z(3) Polyakov loop symmetry [21]. Based on the 
Polyakov loop, given by the product of gauge links in the temporal direction, 

1 N t 

L(n) = -Tr J] ^o(n,n ) , (3.1) 

n =l 

which transforms under Z(3) as L(n) — ► zL(n), one can construct a standard set of observables and 
order parameters. In particular, (L) = in the Z(3)-symmetric, confined phase, and non-zero in the 
spontaneously broken, deconfined phase. 

As explained in the Introduction, one expects this phase transition to persist in the presence of 
sufficiently heavy dynamical quarks. However, these dynamical quarks explicitly break the quenched 
Z(3) symmetry, favoring the real Z% branch and inducing a non-zero Polyakov loop (L) > even 
in the confined phase. The Polyakov loop and associated observables can no longer serve as order 
parameter. The phase transition will be identified as a discontinuity in these observables, from one 
non-zero value to another, in the thermodynamic limit. 

Furthermore, as also explained in the Introduction, this first-order phase transition is expected 
to become weaker and eventually disappear, as the quark mass is reduced. The correlation length 
at criticality will correspondingly increase, and eventually diverge at the end-point of the transition 
line where the transition becomes second-order. For yet lighter quarks, no singularity appears even 
in the thermodynamic limit, and a simple crossover occurs in (L). The distinction between weak 
first-order, second-order, or crossover behaviour requires lattice sizes at least comparable with the 
critical correlation length. However, this requirement can be lessened if one compares results for 
various small size lattices with a finite-size scaling ansatz. This is how the first-order nature of the 
quenched transition was first ascertained Q, and this is how we proceed here. 

Our strategy is to vary /3, for each quark mass, for three spatial lattice sizes and to look for the 
following signals: 

• Coexistence of the two phases: A distinctive feature of a first order transition is phase coexistence 
and on a finite lattice we look for tunneling between the confined and the deconfined phase which 
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is observed over a small temperature range around the critical temperature. As the size of the 
lattice increases tunneling is exponentially suppressed but for the lattice sizes studied here 
enough tunneling events are observed to enable us to study the double peak distribution for the 
norm of the Polyakov loop defined as 

tt = i$>(n) (3.2) 

n 

with V the spatial volume. For a first-order phase transition, the location of the peaks should 
be fixed as the volume increases, while their width should decrease like V~ l . For a second- 
order phase transition or a crossover, the peaks should merge as the volume increases, and the 
tunneling signal should disappear. 

Deconfinement ratio: The deconfinement ratio is defined by 

3 1 , . 

P=- 2 P~- 2 (3-3) 

with p the probability for the complex Polyakov loop to be within 20 degrees of a Z% axis. 
Therefore if the Z{2>) symmetry is unbroken we find p = 1/3 and p = 0, whereas if it is broken in 
such a way that the Polyakov loop is distributed near one axis (the real axis) we have p = 1 and 
p = 1. The value of p = is only obtained in the quenched case where the Z(3) symmetry is 
exact. With dynamical quarks we only look for a discontinuity across the phase transition. On 
a finite lattice, the discontinuity is smoothed out. (p) varies abruptly over a range in (3 which 
shrinks as V , and the slope dp/d(3 is maximum at the critical (3. The intersection of curves 
p(P) for various lattice sizes may also give information on (3 C . [] 

The real part of the Polyakov loop after projection on the closest Z% axis: This observable is also 
an order parameter for the quenched phase transition. It behaves just like the deconfinement 
ratio in the presence of dynamical quarks. 

The peak value of the susceptibility: The susceptibility measures the fluctuations of the Polyakov 
loop. We consider the behaviour of 

XL = vm\ 2 )-mf) (3.4) 

which diverges at criticality for a first-order phase transition in the thermodynamic limit. On 
a finite lattice, instead of this 5 function behaviour, the peak value of \L is proportional to the 
volume V, while the width of the xl distribution scales like V^ 1 and its peak may shift like 



V^ 1 [22]. The scaling behaviour of the susceptibility changes for a second order transition: the 
peak value of xl becomes proportional to V a with a < 1. For a crossover behaviour where, 
even in the thermodynamic limit, there is no discontinuity in the thermodynamic functions, the 
peak value of xl remains constant with the volume. 



One complication here is that (p) in the confined phase depends on the spatial volume V. 



7 



4. Results 



In Table |l| we collect the parameters of our simulations. Adjusting the number of bosonic fields, n, 
so that the acceptance is about 2/3, we can make two observations: Firstly, for a fixed quark mass 
we find that n grows logarithmically with the volume for the same acceptance ratio (i.e. for the same 
relative error in the bosonic approximation of the determinant). This is illustrated in Fig. || which 
shows the k = 0.1 data from Table |l|. Secondly, for a fixed volume, n is approximately inversely 
proportional to the quark mass i.e. 1/n is linear in 1/k. This behaviour is displayed in Fig. |3| where it 
holds for a range of k values. This dependence of n on both the quark mass and the volume confirms 
the expectations of e.g. Ref. |23|| . 

Table 1: 

In the first column we give the k values for the three volumes studied. In the other columns n is the 
number of bosonic fields and acc the average acceptance; Ksw gives (in kilo sweeps) the total number 
of thermalized configurations used in the reweighting procedure 



K 


8 3 x 4 12 3 x 4 16 3 x 4 




n/acc Ksw 


n/acc Ksw 


n/acc Ksw 


0.05 


8/0.78 18 


12/0.74 20 


24/0.83 20 


0.10 


16/0.67 45 


24/0.63 50 


32/0.67 37 


0.12 


24/0.74 55 


32/0.67 40 


40/0.69 12 


0.14 


32/0.77 60 


40/0.70 37 


50/0.67 12 



Our results for the observables which we used to decide on the order of the transition are shown 
in Figs. (|j - |7|). In Fig. || tunneling between the confined and the deconfined phase is clearly observed 
for k = 0.1. A similar behaviour is also found in the case of k = 0.05 whereas for k = 0.12 tunneling is 
no longer observed, on our largest, 16 3 , lattice. The double peak structure of \Q\ for k = 0.10 is seen 
in Fig. [5] where it is fitted to the sum of two gaussian distributions in the complex plane, one centered 
at the origin (confined phase) and one centered at a fitted location along the real axis (deconfined 
phase). The second peak seems to approach the first as the lattice size increases, which would indicate 
a second rather than a first-order transition. A double peak structure was also observed for k = 0.05 
but the distance between the peaks remains the same as the lattice size increased. For k = 0.12 the 
double peak was no longer visible. 

In Fig. I we show the deconfinement ratio p as a function of 0, using reweighting |24]], for k = 
0.05,0.10 and 0.12 as well as for the quenched case (k = 0). The latter is included as a check of 
our methodology and statistical accuracy. Our quenched data has similar statistics to our full QCD 
data. The intersection of the deconfinement ratio curves p(f3) for various lattice sizes occurs very near 
the known value of (3 C [|| shown by an asterisk, and around a value 3/4 for the deconfinement ratio. 
The largest lattice (24 3 ) results, actually obtained by reweighting data from our effective model (see 
Section 5), help determine the critical point more accurately, but are not necessary to ascertain the 
transition itself. A qualitative change is seen to occur around k = 0.10, namely the curves for various 
lattice sizes stop intersecting. In fact for k = 0.14, which is not represented, the deconfinement ratio 
always stays very near 1. 

A similar qualitative change is visible in Fig. 0, which shows the susceptibility \L divided by 
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Figure 3: (a) The number of bosonic fields, n, is plotted vs the logarithm of the volume for k = 0.1. 
(b) 1/n vs 1/n for the 8 3 x 4 lattice. 
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Figure 4: Monte Carlo history of the Polyakov loop for k = 0.10 on lattices 8 3 x 4, (3 = 5.670 (upper), 
12 3 x 4,/3 = 5.670 (middle) and 16 3 x 4,/3 = 5.660 (lower). 
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Figure 5: Double peak structure of |0| for k = 0.10 fitted to a sum of two gaussians, centered at the 
origin and at an adjustable position on the real axis. 

the volume, again reweighted as a function of (3. For k = and 0.05, the peak of xl/V is almost 
independent of V. In other words, xl diverges like V, as befits a first-order transition. For all 
the higher k's, xl/V decreases as the volume grows, indicating that the first-order transition has 
disappeared, or that the asymptotic behaviour has not yet set in. 

We use the peak of the Polyakov loop susceptibility xl to determine the pseudo-critical coupling 
(3 c (k). The dotted lines in Fig. [7| delimitate a band where xl is near its peak value on our largest 
volume. These bands, also shown in Fig. g display the nice agreement between this criterion and the 
crossing of the deconfinement ratio lines when that takes place. We have also looked at the average 
of the real part of the Polyakov loop Q after projection to the closest Z(3) axis, and found that the 
position where data from different lattices cross is again within the same error band. Our values for 
(3 c (k) are collected in Table || and plotted in Fig. ||. Comparing them with available j3 c values for 
Nf = 2 [pfifl , also included in Table we see that the effect of one dynamical quark on the value 
of P c is to shift it from the quenched result of /3 C = 5.6923(4) by about half the amount of the shift 
produced by two degenerate quarks for the same k value. Not surprisingly, when the dynamical quarks 
are heavy, their ordering effect grows linearly with the number of flavours. One could therefore use 
our Nf = 1 data to complete the blank Nf = 2 entries in Table ^. 

The volume dependence of the peak value of xl is more clearly displayed in Fig. |9|a. The lines 
shown are best fits to the form V a . As discussed in Section 3, values a = 1, 0<a<l and 
characterize a first-, second-order transition or crossover respectively, up to finite size corrections. For 
k = 0.05 the best fit yields a = 0.96(4), to be compared to the pure gauge value a = 1.02(4) for which 
the data are plotted in Fig. |^b. The transition is first-order. For n = 0.14 on the other hand we find 
a = 0, a clear signal of crossover behaviour. For k = 0.10 and k = 0.12 the situation is less clear. The 
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k=0.0 /c = 0.05 




5.67 5.68 5.69 5.70 5.67 5.68 5.69 5.7C 




Figure 6: The Deconfinement ratio for n = 0.05, k = 0.10 and k = 0.12 for three lattice sizes. The 
dashed lines give the error band for the critical value as determined from the maximum of the 
susceptibility (see Fig. |7|). The asterisk denotes the quenched result for f3 c from ref. [3]. 



Table 2: 

The pseudo-critical (5 value is given for the four k values studied. The two-flavour results are from 
ref. [25] and are given for comparison. The quenched result from [0] is (3 C = 5.6923(4) 





K 


0.05 


0.10 


0.12 


0.14 


N f = l 


(3c 


5.692(2) 


5.660(3) 


5.630(5) 


5.59(1) 


N f = 2 








5.58(2) 


5.46(2) 
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Figure 7: The susceptibility per unit volume for the quenched theory (k = 0) and for k = 0.05, 
k = 0.10, k = 0.12 and k = 0.14. The asterisk marks (3 C [3] for the quenched theory. 
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small value of a = 0.22(3) at k = 0.12 favors the crossover region. For k = 0.10, a = 0.56(3), and we 
must be quite near the second-order end-point. 
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Figure 8: The pseudo-critical deconfinement coupling [3 C as a function of n Nt (open squares) found 
from the analysis of the QCD data with one dynamical quark species (here Nt = 4). The dashed line 
is to guide the eye. The results for the same quantity obtained for our effective model (Section 5) are 
also plotted as a function of the field strength h. The asterisk shows the second-order end-point. The 
QCD data are shifted upwards by 1QN c k a to take into account an effective change in (3 coming from 
expanding the fermionic determinant up to order k 4 . 

To summarize our findings: 

• For k = 0.05 all our observables (tunneling, intersecting deconfinement ratios p(/3), linear de- 
pendence of the susceptibility peak on V) behave quite similarly to the quenched case and 
consistently point to a first-order transition. 

• For k = 0.12 and 0.14, all our observables (no tunneling, no crossing of the deconfinement ratios, 
growth of the susceptibility peak slow or non-existent) point to a crossover. 

• For k = 0.10, the situation is less clear. On the largest lattice, we do not see the signals of a 
phase transition: tunneling is suppressed as compared to ft = 0.05, the peak in the Polyakov loop 
distribution moves towards the origin, the deconfinement ratio does not intersect the smaller 
volume curves, and the peak of the Polyakov loop susceptibility is not much larger than the 
peak value for the next smaller volume. Based on this largest volume, we consider k = 0.10 as 
being close to the end-point, but rather on the crossover side. 

A more precise statement would require simulations on larger lattices. Larger volumes together 
with a more refined finite size scaling (FSS) ansatz would lead to a precise determination of the end 
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Figure 9: (a) The volume dependence of the peak of the susceptibility for k = 0.05, k = 0.10, k = 0.12 
and k = 0.14. (b) The volume dependence of the peak of the susceptibility in our effective model 
(Section 5) for h = 0.01, h = 0.02, h = 0.025 and h = 0.05. A couple of points are slightly shifted in 
volume for clarity. 

point. We carried out this approach in the context of the effective model discussed in detail in the 
next Section. 



5. Effective model 

It is clear from the above results that locating the end-point (/3 ep , n ep ) accurately requires very large 
lattice sizes, beyond our computer capabilities. We thus adopt a different approach: We consider an 
effective model, in the same universality class as full QCD but cheaper to simulate. By simulating 
larger systems, we locate the end-point in the coupling plane of this model. Then we map this 
end-point back onto full QCD. 

The simplest prototype of our universality class is the three-dimensional three-states Potts model 
in an external field h: 

5 = -V l3Rez*zi-hJ2 Rez n h >® (5-1) 



E 

<nl> 



with z n an element of Z(3). Dimensional reduction at high temperature reduces QCD to this model. 
Quenched QCD maps onto the Z(3)-symmetric h = version, which undergoes a first-order transition 
at j3 c rj 0.55053 [ J26| | . Full QCD reduces to the h > version, where the strength of the Z(3)-breaking 
field h grows with the inverse quark mass. That model has been studied by DeGrand and DeTar 
[15]. They found a first-order transition line in the plane (f3,h), ending at (f3 ep ,h ep ) (see Fig. [l|). 
The second-order end-point critical field h ep was evaluated in the mean-field approximation, yielding 
the small value | log 2 — | ~ 0.018, and by Monte Carlo where the estimate 10~ 3 < h ep < 10~ 2 
was obtained. Unfortunately the mapping from h ep back to a quark mass is qualitative rather than 
quantitative. This is why we have to turn to a more complicated, four-dimensional model. 
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Table 3: 

Non - perturbative determination of h by fitting to the QCD data. The value of h listed here gave 
the best fit to the QCD data of the |0| distribution. 



K 


0.05 


0.10 


0.12 


h 


0.005(4) 


0.026(2) 


0.06(1) 



The starting point is the fermionic determinant which can be expanded into loops yielding 

1 



det(l - kM) = exp{-^^-Tr{M 1 )) (5.2) 

l 



This expansion will converge quickly for the range of k's under consideration. The loops which are 
relevant for the breaking of the Z(3) symmetry are the ones that wind around the time direction. 
Among them the shortest and most important is the Polyakov loop eq.([0]), which carries a coefficient 
jf t (2n) Nt on a lattice of time dimension Nt- Our effective model tries to incorporate higher-order 
effects with an effective Z(3)-breaking field h, and an action 

S cS = S g [U] - h( K ) ^ TrL(x) (5.3) 

X 

Besides the leading hopping parameter expansion above, which should be accurate at small k, the 
mapping h{n) can be obtained in other ways. One way is to resum all spatial hoppings within the 
free approximation, so that 

h(K) = N t jT ^0 Tr [(1 + 70 ) k G(q) ] N * (5.4) 

where G(q)~ 1 is the free spatial quark propagator. This mapping incorporates the divergence of h when 
k — > k c , with k c = 1/8. This result can be further refined by considering tadpole-improvement, such 
that k in (|5.4| ) is substituted by n{plaq) 1 ^ . Finally, a non-perturbative approach consists in finding 
the best match h(n) between the distribution of the magnitude of the Polyakov loop (eq.(p.2j)) in 
full QCD and that in the effective model. The quality of such a matching is illustrated in Fig .[10|- At 
lower k or h, this non-perturbative matching becomes more sensitive to statistical fluctuations in the 
exploration of Z(3) sectors, which can mimic a symmetry-breaking term. In that regime the other, 
analytical methods work better. Our non-perturbative matching results of the whole distribution 
are listed in Table [|. (Matching the first two moments of the distribution gave consistent results). 



Fig. 11 compares the various mappings h(?z) considered here. 

Now equipped with a reasonably accurate correspondence between h and k, we can look for 
the critical line (J3, h) and its end-point in our effective model. We simulate lattices of spatial size 
8 3 ,12 3 ,16 3 and 24 3 with the temporal direction fixed at Nt = 4 as in the case of the full QCD 
simulations. We obtain data at various values of h = 0.01,0.02,0.025 and h = 0.05 with typically 
10000 — 20000 thermalized configurations. 

First, we perform the same analysis, using the same observables as in Section 4 for full QCD. At 
each value of h, we identify the pseudo-critical coupling f3 c (h) and obtain the points included in Fig. ||, 
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Figure 10: Best fit of |fi| from full QCD for k = 
0.1 (dashed line) to the data obtained from the 
effective model using reweighting (solid line for 
h=0.026). 
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Figure 11: The strength of the Z{3) breaking 
term versus K 



which for small k follows closely the corresponding curve obtained in full QCD. In Fig. (gb) we show 
the dependence of the peak of the Polyakov loop susceptibility xl on the spatial volume for various 
values of /i(k). The behaviour is again very similar to that of full QCD, indicating a weakening of the 
first-order phase transition as h increases, and probably a crossover behaviour for h > 0.03. 

Now we would like to determine as precisely as possible the end-point (f3 ep ,h e p), and check the 



scaling behaviour in its vicinity. In ref. [27] the same issues were addressed in the context of the 
electroweak theory, and a reliable numerical procedure was presented which yielded impressively 
accurate answers. The procedure we follow is very similar, although not quite identical. From our 
Monte Carlo data we obtain the joint probability distribution of the plaquette and the real part of 
the Polyakov loop (shown in fig. ^). The principal axes of this distribution, which diagonalize the 
correlation matrix of the two observables, are identified as the magnetization-like variable M and 



the energy-like variable E. (The rotated distribution displayed in fig. 13 shows the double peak 
distribution in the M - like direction.) After subtracting the averages from the new variables M and 
E such that (M) = (E) = 0, rescaling them such that (M 2 ) = (E 2 ) = 1, and reweighting in /3 and 
h, we obtain the probability density P(M, E) shown in Fig. |14| at the end-point. We clearly see that 
the marginal distribution P{E) will show a single peak, while P(M) will have a double peak. In an 
infinite volume the distribution P{M) would be symmetric. Comparing the P(M, E) distribution of 



Fig. [T| with the 3d results of Fig. 3 Ref.|2j for the Ising model and the 0(2) and 0(4) models, our 
result is in closest agreement with the P(M, E) distribution of the Ising model. 

To determine the end-point we need to perform a finite-size scaling analysis. For each lattice size, 
we reweight in [3 and h until we minimize the asymmetry of the M distribution, by requiring the 
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vanishing of the third cumulant: 

WW = (5 - 5) 

This determines a line in the plane ((3,h), which must go through the infinite- volume critical point 
((3 ep , h ep ), up to statistical errors. The intersection of these lines obtained for various lattice sizes 
therefore determines the critical point (/3 ep , h ep ). Parametrizing these lines by the value of h, we show 



in Fig. |15| the variation of the fourth magnetic cumulant 

(M 4 ) 
(M 2 ) 2 



(5.6) 



along these lines, for the three largest lattices considered for Nt = 4 as well as for three lattices for 
Nt = 2. The simulation of the effective model at Nt = 2 using spatial sizes 8, 12 and 16 can be used 
as a check of our procedure in a situation where significant deviations from scaling are not expected. 
The three lines for N = 2 cross at a point which gives the critical end-point. The fact that the three 
lines meet at a point (within statistical errors) indicates stability against deviations from scaling. 
Note that we made no assumption about the universality class in our determination of the critical 
point. Because of its generality, and because of its relative statistical robustness in our case, we favor 
the method we just described to find the end-point over the method used in Ref. [p7f| , which fixes the 
fourth cumulant eq.( |5.6| ) to its Ising value. The value of ((AM) 4 ) / ((AM) 2 } 2 at the critical end-point 
can now be used to check the universality class of the effective model and thus the universality class 
of QCD. From Fig. |l|(a) at h ep = 0.06 we find ((AM) 4 ) / ((AM) 2 ) 2 ~ 1.67(5) in good agreement 
with the corresponding value of 1.604(1) for the Ising modelQ. The value of /3 ep determined from the 
crossing of the lines parametrized with the value of (3 is /3 ep = 5.047(1). For Nt = 4 the lines for 
our three largest volume converge at h ~ 0.009. Scaling violations are still seen for the lattice with 
spatial extent L = 12. Thus here the end - point is determined less accurately as compared to the 
Nt = 2 case, although we have increased our statistics five fold using ~ 100, 000 configurations for 
the two smaller lattices and ~ 50, 000 for the largest lattice. Since the point of intersection for the 
two largest lattices is within statistical errors consistent with h = 0.009 we take h ep ~ 0.009(1) as the 
best determination of the critical point from these data. Having determined the critical point we can 
look at critical exponents. The magnetic susceptibility at h ep should scale like 



XM = V((AM) 2 ) oc L" 3+ ^ (5.7) 



We plot the \M distributions in Fig. 16 for our three lattices after scaling with L~ ' 8 . They are seen 
to nicely fall on top of each other. This scaling thus yields a rough estimate of the exponent j/v ~ 2.2 
to be compared with a value of 1.96 for the Ising and 0(2), 0(4) model. To determine exponents like 
a/v which can pin down the Ising model universality class we need to consider the scaling behaviour 
of the E - like variable, which is smaller that the M - like variable by at least two orders of magnitude. 
Scaling based on the E - like variable was found to be unstable within the statistical precision of our 
present simulation. 



Finally, we comment on the Polyakov distribution displayed in Fig. 14. It is clear from this figure 
that even at the end-point h ep = 0.009, the Polyakov loop distribution still shows a marked double 

2 For the 0(2) and 0(4) models ((AM ) 4 }/((AAf ) 2 ) 2 = 1.233(6) §§] and 1.092(3) (2§ respectively i.e. clearly lower 
than our value of 1.67(5). 
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Figure 12: (S cS - (S e ff)) vs (S g - (S g )) for 5000 
configurations on 24 3 at h = 0.01, j3 = 5.680 



Figure 13: (S cS - (S cS )) vs (S g - (S g )) for 5000 
configurations on 24 3 at h = 0.01, f3 = 5.680 
after a rotation 



peak structure. In ref. [|(| where the same effective model was studied, the criterion used to identify 
the end-point was to look for the value of h where the double peak structure of the real part of the 
Polyakov loop was no longer visible. This criterion thus led to a much bigger value ~ 0.08 for the 
critical field strength than we are finding here. 
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Figure 14: Normalized probability distribution at the critical point h ep = 0.009 for 24 3 
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(a) (b) 
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Figure 15: The fourth magnetic cumulant versus the field strength h. (a) for N t = 2, (b) for N t = 4. 
The point of intersection of the lines determines the critical point h ep . 




Figure 16: Scaling of the M - like distribution at the critical point h ep and (3 ep for the three largest 
lattices. We have applied a scaling factor of L -0 - 8 to all the data. 
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6. Discussion 



From the previous section we conclude that the end-point of the first-order transition occurs at 
~ 0.009 which maps to k ~ 0.08. The qualitative picture expected for (2+1) flavours in the 
continuum is that there is a finite region in the mass plane {m u ^ m s ) of the two degenerate flavours 
versus the third near the quenched limit where the deconfinement phase transition is first order. Our 
Nf = 1 result corresponds to the infinite mass limit for m u ^. In this limit, which up to now received 
little attention, we find that the first order pure gauge phase transition persists up to k ~ 0.08. Since 
we expect dynamical quark effects to be twice as strong for two flavours (see Table |2|) we can estimate 
the boundary of the first-order region as 

2h eff (k u4 ) + h eff (k s ) « 0.009 (6.8) 



with h e ff{ti) as per Fig. 11 



It is of course of crucial importance to establish the scaling behaviour of our results, and to express 
the end-point parameter as a quark mass, not a hopping parameter. One might wonder if the rather 
small value of K ep we found may not be an artifact of the relatively coarse lattice discretization, 
and if one would recover K ep = in the continuum limit. This would mean that the deconfinement 
transition disappears as soon as dynamical quarks of any mass are allowed. Such a scenario goes 
against the robustness of a first-order phase transition. The first-order deconfinement transition must 
be robust against small variations of the parameters, before ending at a second-order point. We made 
this argument in the introduction to justify the expected phase diagram Fig.[l] in our discrete theory. 
The same argument applies to the continuum QCD theory. Therefore, we expect the deconfinement 
transition to end at a finite quark mass m q . Determining this mass would require a full-blown scaling 
study of 1-flavour QCD, which is beyond the scope of this paper. However, we have considered what 
happens to our effective model as the lattice spacing a changes. We have simulated our effective model 
for Nf = 2 and found h ep ~ 0.06, with scaling exponents consistent with those at iVj = 4. Changing 
Nt = 4 to Nt = 2 amounts to doubling the lattice spacing a. Therefore, the variation of the end-point 
value h ep indicates a scaling behaviour 

h ep (2a)/h ep (a) « 0.06/0.009 oc a 2 - 7 (6.9) 

This strongly supports the existence of a continuum limit for the effective model eq. ( |5.3| ) , with action 

S eff = S[A] -hi d 3 xL(x) (6.10) 

where h = ha 3 in the discretized theory, having a critical end-point at some positive h ep . 

If one believes in the physical nature of the critical coupling h ep , then the corresponding quark 
mass cannot be infinite. For instance, the leading hopping parameter expansion eq. ( |5.2[ ) at constant 
physical temperature T = (iVja) -1 gives 

h{K) ~ AN^ln)^' 1 (6.11) 

If h{n) = h ep a 5 , then this relation is not consistent with lim a ^QK = 0. One is again lead to expect 
persistence of the first order transition, for some range of quark mass, in the continuum one-flavour 
QCD theory. 
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We can get an estimate of the critical quark mass in physical units by using the naive prescription 
m q a = g(i — — ). Using the quenched data of ref. for the pion mass we find k c = 0.1694(2), at j3 = 
5.7 in the zero-flavour theory. On the other hand the SESAM collaboration finds k c = 0.1585(1) |32j 
at P = 5.6 in the two-flavour theory. Going to our value of (3 ep = 5.683(3) will decrease the critical 
value from SESAM slightly. If we neglect this change and take the average between the zero- and 
two-flavour cases, we end up with k c ~ 0.164. Taking the end-point value K ep ~ 0.08 we find a quark 
mass m q a = 3.2 in units of the lattice spacing a. This mass is in fact of the same order as that of the 
doubler modes, which leads to an overestimation of the mass at the end point. Since m q a > 1, finite 
a corrections are expected to be rather large. At tree level the corrected mass is given by 

m„a 

m q a ~ g (6.12) 
y/1 + m q a 

which reduces our estimate by about a factor of two to m q a ~ 1.56. To convert this to physical units 
we use (4a) -1 ~ 220MeV from the deconfinement temperature. This gives m q ~ lAGeV for the 
end-point. 

A critical quark mass of the order of a GeV is in line with phenomenological expectations. The 
pure gauge deconfinement transition is fairly weak, with a critical correlation length £ c ~ (7 — 10)a at 
(5 = 6 pq] i.e. ~ llGeV' 1 or 0(5 a' 1 / 2 ). This is the minimum system size necessary to observe the 
first-order nature of the deconfinement transition. Dynamical quarks introduce a new length scale r c , 
namely the distance where the string breaks, r c ~ 0{2m q / a). Confinement can only be observed up to 
this distance. For very heavy quarks, r c > £ c : string breaking occurs for a large separation, larger than 
the critical correlation length. The passage from confinement to deconfinement will allow "liberated" 
quarks to appear even if their separation is less than r c . This qualitative change signals a phase 
transition. As the quark mass is lowered, the string-breaking scale r c decreases. When r c ~ £ c , the 
passage from confinement to deconfinement does not liberate quarks that were not already liberated 
by string breaking. There is no qualitative change from one regime to the other, and one cannot really 
tell if the system is confined or deconfined: the phase transition has disappeared and been replaced 
by a crossover. This occurs for 

m q ~ 0(5y/a/2) i.e. m q ~ 0(l)GeV (6.13) 



7. Summary and Conclusions 

The multiboson method can be used to simulate an odd number of flavours as well as an even number. 
In this work we have shown that the multibosonic algorithm is well suited for the study of one - 
flavour QCD for moderately heavy Wilson quarks. Using this algorithm we were able to carry out a 
detailed Finite Size Scaling (FSS) analysis to determine the pseudocritical f3(n) line for k values up to 
k = 0.14. We demonstrated that the first order phase transition seen in the quenched theory persists 
when one includes dynamical quarks. Using FSS we showed that the deconfinement phase transition 
gets weaker as the dynamical quark mass increases, then turns into a crossover. In general we find 
that the dynamical quark effects on the phase transition are approximately half those for two flavours. 
Using lattices of sizes 8 3 , 12 3 and 16 3 x 4 within full QCD, we found an end - point around k = 0.1. 
A more accurate determination would require simulation of larger spatial volumes. In order to carry 
out a more refined FSS analysis with larger lattices we considered an effective model in the same 
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universality class as full QCD. In this effective model the effects of the fermionic determinant were 
simulated by an effective Z(3) - breaking field coupled to the Polyakov loop. We studied the phase 
transition as a function of this field strength h on lattices with spatial volumes 8 3 , 12 3 , 16 3 and 24 3 
and found a first order transition that gets weaker as the field strength increases, exactly like in QCD 
as k increases. Performing a FSS analysis using the joint probability distribution of the plaquette 
and the real part of the Polyakov loop we were able to determine, without any assumptions about 
the universality class of the model, the end point of the first order transition line. Matching of the 
Polyakov loop histograms or its first two moments in the effective model to those in full QCD enabled 
us to determine non - perturbatively the correspondence between h and k. In this way we obtained 
the value of n ep ~ 0.08 for the end - point of the QCD first order transition line. Futhermore, the 
FSS analysis of the effective model at the critical point yielded results consistent with our effective 
model (and one-flavour QCD) being in the same universality class as the Ising model. 

We have also studied our effective model on a coarser lattice at Nt = 2, and have found good 
scaling behaviour of the end-point field strength, with similar Ising-like exponents at the end-point. 
This scaling test strongly supports the existence of a critical end-point at a non-vanishing positive 
field strength in the continuum limit of our effective model. In turn, this indicates a non-vanishing 
quark mass for the end-point of the deconfinement transition line T c (m q ) in continuum one-flavour 
QCD. Our findings are in agreement with the idea that a first-order transition is robust against small 
variations of the parameters, before ending at a second-order point. Converting the value of K ep = 0.08 
to physical units can only be done approximately, using quenched and two-flavour determinations of 
k c at similar (3 values. After correcting for finite lattice spacing errors to tree level we obtain an 
estimate of the quark mass at the end - point of m q ~ 1.4 GeV, consistent with phenomenological 
expectations. A simulation of one-flavour QCD closer to the continuum limit would allow tighter 
control over discretization errors, but is beyond our computer resources. 
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